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Viscous corrections to relativistic hydrodynamics, which are usually formulated for small velocity 
gradients, have recently been extended from Navier-Stokes formulations to a class of treatments 
based on Israel-Stewart equations. Israel-Stewart treatments, which treat the spatial components of 
the stress-energy tensor nj as dynamical objects, introduce new parameters, such as the relaxation 
times describing non-equilibrium behavior of the elements Tij . By considering linear response theory 
and entropy constraints, we show how the additional parameters are related to fluctuations of Tij. 
Furthermore, the Israel-Stewart parameters are analyzed for their ability to provide stable and 
physical solutions for sound waves. Finally, it is shown how these parameters, which are naturally 

■ described by correlation functions in real time, might be constrained by lattice calculations, which 
' are based on path-integral formulations in imaginary time. 
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CsJ I The liquid-like quark matter observed at the Relativistic Heavy Ion Collider (RHIC) [l| has inspired a re- 
naissance in hydrodynamic modeling at high energy. However, even though the matter is in a liquid state, the 
extremely short characteristic time and distance scales necessitate the inclusion of viscous effects in hydrody- 
namic models. Viscous effects, or equivalently the effects of non-equilibrium values of the spatial components 
of the stress-energy tensor , are especially important during the first few fm/c of the collision which is char- 
acterized by large velocity gradients and enormous shears. Bulk viscosity might play an important role during 
the hadronization stage, arising from the matter's inability to maintain equilibrium when traversing the phase 
transition [1, i, i, i, Q . 

The standard means for implementing viscous effects into hydrodynamics is the Navier-Stokes (NS) equations 
" ' I 0. However, during the last few years, and especially for relativistic heavy ion physics, increasing attention 

■ has been focused on the Israel-Stewart (IS) equations [1, H, [13, [HI- Whereas the spatial components of the 
I ■ stress-energy tensor are determined by the energy density and velocity gradients in NS theory, in IS treatments 

Qs^ , the components are treated as dynamic objects which decay exponentially towards the NS values. In this 
Cf^ ' paper the IS equations are analyzed in detail, the differences between IS and NS descriptions are explored, 
' , and the degree to which the additional parameters required for IS theory might be microscopically calculable 
I ■ from approaches like lattice gauge theory is investigated. After reviewing the basic motivations and premises of 
f — I IS-based equations in this section, the formulations are analyzed from the perspectives of entropy growth and 
. linear response theory in Sec. |TT1 These perspectives constrain both the form and parameters of IS theory. In 
^ ' particular, the relaxation times, along with the viscosities, are shown to determine both the functional form 
for the relaxation and the statistical fluctuation of the elements of the stress-energy tensor. Some of these 
constraints are also derived from the perspective of linear response theory and are presented in the subsequent 
section. Section IIVI presents a discussion of how one might alter IS equations to be more physically relevant for 
large deviations from equilibrium, and Sec. |V] provides an analysis of the range of IS parameters that provide 
stable and physical sound waves. The prospects for determining IS parameters from lattice calculations are 
discussed in |Vl] . 

The traditional means for incorporating viscosity is through the Navier-Stokes equation, in which the stress- 
energy tensor is modified by the velocity gradients in a linear fashion, 

= P5,, - Tj [d^v, + djv., - (2/3)V • v] - C%V • V. (1) 

Here, it has been assumed that the stress-energy tensor is expressed in the rest-frame of the fluid, where 
V = 0. The pressure P, the shear viscosity 77, and the bulk viscosity C, are all functions of the energy and particle 
densities, though we will ignore the particle-density dependence throughout the remainder of the paper. Since 
the stress-energy tensor is symmetric, the deviation of Tij from the equilibrium value P5ij can be expressed in 



terms of six independent numbers: 
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Here, b refers to the deviation from equilibrium of the trace of Tij and is related to the bulk viscosity, while the 
five values ai , which are related to the shear viscosity, describe the angular anisotropy. In addition to V • t; one 
can also define five symmetrized "velocity gradients" , 



uJi = - {2dxVx - 2dyVy) , 
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[dxVx + dyVy - 2dzVz) , 



^^3 = {dxVy + dyVx), UJ4 = {d^v^ + d^Vx) , uj^ = (dyV^ + d^Vy) . 

With these definitions, the work (in the local rest frame) due to expanding a volume V for a time increment dt 
becomes. 



■^dW — dt ^ TijdiVj = dt 



P(V ■v) + bV -v + ^ai, 



(4) 



and the Navier-Stokes equations assume a simple form, 

6(NS) ^ _^vz;, (5) 

(NS) 

a] = -77WJ. 

Israel-Stewart approaches Q differ from NS in that and b are treated as dynamical objects. Originally, 
the motivation for IS approaches was related to their numerical stability compared to NS equations [9|. This 
property follows from the fact that the equations are hyperbolic, as opposed to NS equations which are parabolic, 
similar to the diffusion equation. Since the diffusion current can move arbitrarily quickly for arbitrarily large 
gradients, the equations were effectively non-causal. Though the violation of causality was never shown to be 
important in the context of RHIC collisions, IS solutions have been considered safer than NS solutions. An 
analysis of the stability of sonic modes of a given wavelength is presented in Section |Vl 

The basic premise of IS solutions is that the offsets of t,j decay toward their NS values exponentially. 
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where a and /? are functions of the energy density. They are usually set to unity, and if the energy density is 
fixed, the forms for a and /3 are irrelevant. Providing an energy dependence to a and P change the equations 
of motion for large velocity gradients. For instance, with a large velocity gradient a system will expand by a 
given fraction in a time inversely proportional to the velocity gradients. If the expansion time is much less than 
the relaxation times and T{,, the ratios a.^/a and b/p will be effectively frozen for very rapid expansions. As 
an example, if a were set equal to P, a very rapid expansion would freeze the ratios Qi/P, whereas setting 
a = I would lead to a freezing out of the variation from equilibrium a^. For a rapidly cooling system where P is 
rapidly falling, the choice of a = P might keep the components of Tij positive, whereas the choice a = 1 might 
lead to negative elements. 

As an aside, we mention that the derivatives D /Dt in Eq.s © include the effects of rotational flow 10]. For 



a tensor r, 
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This differs from the usual notation where D/Dt also includes the term v • V. However, that term is explicitly 
zero in the rest frame of the fluid. If one were solving the equations of motion for f2 ^ 0, one would need to 
rotate the stress energy tensor at every time step consistently with Eq. ([7]) . This would mix the five elements , 
but would leave b unchanged. Since is proportional to velocity gradients, and since is also proportional to 
viscosity gradients in NS theory, corrections for rotational flow are considered a "second-order" correction from 
the NS perspective. Thus, one expects rotational effects to be negligible except for those instances where there 
exists both rotation and significant values for the anisotropics a^. Henceforth, we will neglect rotation, though 
we emphasize that it might become no n- negligible for non-central collisions, or for fiow away from mid-rapidity. 

The parameters of IS hydrodynamics are the viscosities, ry and B, the relaxation times, Ta and Tb, and the 
scaling functions for the exponential decays, a and f3. Each of these six parameters can be functions of the energy 
density. However, as will be seen in the next section, the constraint that the entropy cannot fall constrains the 
number of parameters to four, and relates them to the statistical fiuctuations of Tij. 



II. ENTROPY CONSTRAINTS 



Entropy must rise, and since the IS equations describe entropy production, the IS parameters are constrained. 
Furthermore, the entropy should depend not only on the energy density, but also on the deviation of Tij from 
equilibrium. In this section we will show how the entropy constraint reduces the number of IS parameter from 
six to four, and relates them to the statistical fluctuation of r^ , which in turn is related to how the entropy is 
affected by non-zero values of and b. 

For a fixed energy density, should equal P if the entropy is maximized. For small variations, and 6, the 
entropy penalty will grow quadratically with a and b. 
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where Ca and ab are functions of the energy density. Since the probability for a given fluctuation of the stress 



energy tensor occurs with probability e*- 



, the variances for ai and 6 in a volume V are 



V 



(9) 



The fact that the fluctuations of ai and b fall as 1/yV is expected since the correlations persist over a finite 
domain and the number of independent domains should increase linearly with V . Equivalently, one would expect 
the fiuctuation of the quantities aiV and bV to increase with ^/V . 
The change in entropy is then 



dS = dS, 
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The change in the equilibrium entropy for an expanding volume element is 

d^cquii = (dE + PdV)/T. (11) 

where P and T are functions of e assuming equilibrium. 

After using the fact that dV = Vdt{W ■ v), and using Eq. ([4]) for the differential work, which equals dE, 
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As expected, entropy production vanishes for = 6 = 0. 
Inserting the IS equations of motion in Eq. ([6]) into Eq. 
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For the entropy to rise regardless of the values of a and b or the velocity gradients, the terms linearly 
proportional to ai and b must disappear as well as those linearly proportional to the velocity gradients and 
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V • V. Furthermore, since a and /3 are functions of the energy and can change arbitarily quickly for arbitrarily 
rapid expansions, ensuring that the entropy always rises requires stating that both a/aa and /3/crf, are constants, 
which can be set to unity without changing any behavior, 

a = (Ta, /3 = o-fc. (14) 

If the terms linear in a and b are to never contribute to the entropy regardless of the values of a and b, 

cr„ = , CTb = — . (15) 

Thus, the viscosities and relaxation times uniquely determine the scaling functions for the exponential decays, a 
and /3, which are equivalent to the fluctuations of CTq and df, at fixed energy. This latter expression was derived 
from the perspective of the Boltzmann equation [9J. This derivation is both more general, and provides the 
equivalence between the fluctuations and the scaling functions. 

After the parameter constraints are enforced, the growth rate for the entropy has a simple form, 

V dt C ^ r] 

After enforcing the NS equations in Eq. ([5]), this reproduces the usual textbook result that the rate of entropy 
production is proportional to (V • v)^ and uif. The crucial difference between the entropy production rate in 
Eq. and the equivalent NS result becomes important in the limit of a very rapid expansion. The time 

required for the dimensions of the volume to increase by a specific amount will be proportional to 1/(V • v) or 
l/uJi- Considering a time step, 

At = 7/V • V, (17) 

a fluid element will increase by a factor, 

Vit + At) ^Vit)e'''. (18) 

For a small 7, the equations of motion for ai/aa and b/at in Eq. © show that in the limit of large velocity 
gradients they change by an amount: 

A{b/ab) = 7C/T6, Aajaa = ,^'^\ , (19) 

(V • v)Ta 

which, as Eq. (|16p shows, results in the entropy production proportional to 7. Thus, in the limit of infinitely 
fast velocity gradients, the net entropy production approaches a fixed value depending on 7. In contrast, the NS 
production rate is proportional to the velocity gradients squared, so that the net entropy produced for a fixed 
expansion factor 7 is divergent as V • d — > 00. For the opposite extreme, of a small velocity gradient, the values 
of b and in Eq. (fTB|) will relax to the NS values, and both the IS and the NS results for entropy production 
will be identical. 



III. FROM THE PERSPECTIVE OF LINEAR RESPONSE THEORY 

According to hnear response theory [1, [13, the classical hmit of the Kubo formula gives. 



V=^l dro I d\{n,{0)n,ir)), i^j (20) 
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dro J d^r {ai{0)ai{r)) , any i. 



= J^l^dro J d'r mb{r)) . 

Here, the averages involve summing over states at fixed energy. If the energy is allowed to vary, as in a grand 
canonical ensemble, one must subtract the contribution from fluctuating energy, c^{SE'^), from the expression 
for the bulk viscosity. 
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For a fluctuating field that loses correlation exponentially, ~ e 



■n = (af) , ai = ^ J (fr a,{r), 



Given that the entropy penalty for small is, 



(fr 



one can identify the variances of a and 6, 



(21) 



(22) 



The 1/V factor was expected given that there are no long-range correlations in Tij. 
Plugging this into the expression for the viscosities, 



T ' ^ T 



(23) 



(24) 



which matches the result of the last section. 



IV. LARGE DEVIATIONS FROM EQUILIBRIUM 

As shown in the previous sections, IS formalism handles arbitrarily large velocity gradients, but remains an 
expansion in the deviations and h. The non-equilibrium values a/ua and b/ab in Eq. ^ decay towards 
values which can be arbitrarily large for arbitrarily large velocity gradients. However, this become unphysical 
in some cases as the elements Tij become arbitrarily large, or strongly negative. Depending on the source of the 
viscosity, it might be more physical to impose a limit on the deviation of from equilibrium. For instance, 
if the shear viscosity is caused by a finite mean free path, the kinetic pressure should remain positive, whereas 
if the dynamics are those of classical electric fields, the stress-energy tensor should be confined to the region, 
— e < Tij < e. Here, we provide a simple recipe to enforce such constraints, within the IS picture, that do not 
violate entropy constraints. 

First we consider the case of the bulk viscosity, and assume that there exist some physical constraints that 
enforce |5| < 6max- To dynamically limit 6, we alter the equations of motion, 

^ = -[.x-(CV.z;)M]/r,, (25) 



dt 



b = 6maxtanh 



For small b these become identical to the previous IS equations of motion. The equations for entropy production 
become. 

Again, the term proportional to V • u must vanish if the entropy is to always grow despite the sign of the velocity 
gradient or the value of x. This determines ds/dx, 

^=-b'^---, (27) 

dX 1 Th (Tfc 

where Eq. (|15p was used to simplify the prefactor. Integrating to find the entropy density, 

s = seq - %^ In [cosh(abx/6^ax)] . (28) 
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The rate of entropy production is then manifestly positive, 



= T^^x = xT^^ tanh(cr6x/6niax), (29) 
V at atn crtn 



which by virtue of Eq. (p4)) becomes 

T dS bxat xab bb^^^ \^-1(uik \ r-in\ 

V~dt ^ ^ ~ '^'^^^Z — ~ (o/Omax) ■ (30) 

Again, for small b this gives the identical behavior as the previous sections. 

A similar procedure can be followed for the shear terms to enforce the constraint, 



by assuming equations of motion, 



^ = -—iyi-V^i/ca), (32) 

at Ta 



Omax tanh I ) , y ^ \l yl + vl, 



Qi — a—, a = \l -\- a%. 
y V 

After enforcing the constraint that the entropy grows regardless of the value of uji , one comes to the result for 
both shear and bulk deviations from equilibrium. 



rd5 ^ ^^j^^]^(-2;(T6/6niax) + ^^tanh(?;a-a/amax) 

V at C, (Ta 



(33) 
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= Seq - In [cosh((7f,a;/5niax)] - In [cosh(cray/amax)] ■ 



The expression for entropy generation is equivalent to that described in Eq. (fT6|) in the limit of small / flmax 
and 6/6max- 

This procedure is more phenomenology than theory. Whereas the coefficients ry, C, aa and Ub can all be 
formally expressed in terms of microscopic correlations, the additional parameters flmax and femax, as well as 
the formalism itself, are chosen at the discretion of the modeler. In fact, one might determine the parameters, 
or alter the formalism, based on a numerical comparison of the IS evolution described above with a particular 
microscopic model. For instance, one might run a model based on particular color-glass assumptions in a 
simplified geometry, then tune the IS approach to match the behavior of the stress-energy tensor. The IS 
model could then be considered a surrogate for the color-glass model, but would more naturally couple to the 
more theoretically justified hydrodynamic treatment of the later stage, and also more naturally incorporate 
complicated three-dimensional evolution. 



V. STABILITY OF ISRAEL-STEWART SOLUTIONS 

One of the often-stated advantages of IS equations derives from the numerical stability of the approach . 
Furthermore, the solutions to NS equations can contain super-luminar transport. Here, solutions to IS equations 
are presented for small perturbations. After solving for the angular frequency, uj{k), the group velocity dio/dk 
is plotted. In the limit of the relaxation times becoming zero, the IS and NS solutions become identical, and 
high momentum modes become super-luminar or unstable. For longer relaxation times, the solutions are well 
behaved, and a characteristic scale for the transition from poorly- to well-behaved can be expressed in terms of 
the viscosity, speed of sound, and enthalpy. 

To solve for the sound waves, we first consider perturbations of the energy density 6e, the velocity v and the 
deviation of the stress-energy tensor, Stzz ■ 

Se = Ee"^'+'''', (34) 
6tzz = Ae''^'+''''- , 
Vz = Ve''^*^'''' 
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The IS equations of motion for small E, A and V are: 

^'if--iiclSe + Sr..), (35) 

— = (r,. + C— 

dt Ti, \ dz 

dSe dv 

It " ' dl' 

Using the plane- wave form of Eq. p4p . and assuming = 0, the IS equations of motion become: 

ujE = ~hkV, (36) 

n Tb 

kcl ^ k , 
ojV = — T^E- -A. 

h h 

One can then solve these for x = iuj/kc and obtain a cubic equation, 

+ Cx^ + {l + -f)x + C = 0, (37) 



C EE 

7 = 



1 



kcsU ' 
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These approach the NS equations when — > 0. 

The solutions to Eq.s p7|) are shown in Fig. [TJ The three panels illustrate the solutions as a function of the 
wave number k for the real and imaginary parts of lo, and the group velocity duj / dk. To provide a more general 
expression, frequencies are given in units of the inverse characteristic time, 

'^'^har = (38) 

which has dimensions of time provided Cs is given as a fraction of the speed of light. This time characterizes 
the strength of the viscous term. The three curves then depend only on the ratio of the relaxation times to 
the characteristic times, T{,/Tchar- For small relaxation times, Lij{k) becomes negative and diverges, leading to 
arbitrarily large negative phase velocities. For r^/Tchar >^ 0.16, the divergence disappears, but solutions remain 
rather peculiar in that the speed of sound remains super-luminar for certain wavelengths. For larger relaxation 
times the behavior seems normal, and the phase velocities only modestly differ from the nominal speed of sound. 
Not unexpectedly, in the same region where the phase velocities most greatly differ from the speed of sound, 
the imaginary part of cj, which describes the decay of the sound wave, becomes large. 

Since the solutions uj[k) were generated by solving a cubic equation, there arc three solutions. Two of them 
correspond to plane waves moving in the positive and negative directions, and are described by the solutions 
plotted in Fig. [TJ For the third solution, the real part of lo is zero, and the waves do not propagate, and only 
decay. 

The solutions presented here are for the case with only bulk viscosity. One can solve the analogous set of 
equations as ((36)l for having shear viscosity, but no bulk viscosity. In that case, Eq.s ([2]) give 6tzz = (27-^^)02, 
and given the Eq.s ([3]), the IS equations of motion for 02 are 



da2 _ 1 

dt Ta 



a2 + {2/Vi)iidzv} ■ (39) 



Finding uj{k) for the sonic waves gives the same result as shown in Fig. [TJ only with the definition, 

477 

re.. = (40) 

If one were to simultaneously consider both shear and bulk terms, one would have to simultaneously solve four 
equations of motion, rather than the three given in Eq. (|36p . This would lead to a quartic equation, unless the 
two relaxation times were equal. 
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FIG. 1: (color online) As a function of the wave number k, the real and imaginary solutions for u}{k) = E-[-iT are shown 
in the upper two panels, while the phase velocity, dw/dk, is shown in the lower panel. By scaling the wave number by 
the characteristic time, defined in Eq. (|38|l . the solutions depend only on the ratio of the relaxation time rt/rchai- For 
Tt — > (solid line), the solutions approach Navier-Stokes solutions, which give peculiar behavior for wave numbers larger 
than l/csTchar. The behavior ameliorates for higher relaxation times, Tb/rchar ~ 0.2 (dashes), Tb/rchar ~ 0.4 (dotted), 
Tb/Tchar = 1.0 (dot-dashed) . 



VI. DETERMINING COEFFICIENTS FROM PATH INTEGRALS 

The transport coefRcients used in IS theory can all be expressed in the terms of correlations. Using the full 
quantum expressions [3], 



c 



^ / d\ {5{T,.M + Tyy{0) + T,M)5T,x{t^0,v)) 



(41) 
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d^x xo {[S (r,,(0) + Tyy{0) + T,,(0)) , ST^^ix)]) 



where the brackets describe a thermal average, (A) = TiAe~^^ , the time dependence of the operators are 
defined by A{t) = e*^*^e~*^*, and the deviations from vacuum expectations are denoted by SA = A — {A). 

CoefRcients can be calculated reliably in perturbation theory at high temperature when the coupling is small 
[13, Unfortunately, lattice calculations, which should be valid for large coupling, do not yet provide reliable 
values for the viscosity The difficulty with calculating the correlations above from lattice calculations 
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derives from the fact that the operators are evaluated for real times, whereas the path integral used in a thermal 
trace samples a path over imaginary time, 

Z = J d$ exp^^-i J d'^xjC{^,d^)^ , < xo < i(3 (42) 



3x]y3z3l3 



Herre, the lattice points are denoted by jx,jy,jz,ji3- In order to represent a trace, the fields are cyclic in the 
time variable, 

<^(0,f) = <^(i/3,f), (43) 

while each of the spatial dimensions is confined to a length L. The expression becomes exact in the limit that 
the lattice points j become infinitely numerous. Although the expression above assumes scalar fields, more 
complicated expression for Fermi or for vector fields will not change any of following derivations. 

Correlations, similar in appearance to those used to calculate the viscosities, can be expressed in terms of 
path integrals, 

V{A{0)B{t)) = ^ J d(f A{0)B{t)expS^-i J d'^x£{<p,d(l))^ , (44) 

where the operators A and B can be expressed in terms of the fields, assuming that t corresponds to a point 
on the path, and tht P orders A and B according to the ordering on the path. For instance, A might refer to 
J d^xTxy, where the stress-energy tensor can be expressed in terms of the fields, 

T^f'ix) = -£{x)g^^ + —dpcl>. (45) 

The notational burden can be reduced by referring to correlations of an operator with itself as 

Gr^it) ^ {A,{Q)A^{t)), (46) 
A = ^ / d\ T,y{r), 



TxX + Tyy + Tg^ 2jrp X 

o Cg(Joo - e) 



In the last term, the subtraction of c^(T'oo — e) enforces the correlator to vanish for large times. Otherwise, 
energy conservation would result in the fluctuation of the energy density (which is related to the specific heat) 
causing a correlation between the elements at infinte relative time. With these definitions the transport 
coefficients become, 

= G,(0), (47) 

~i r°° 

= Y y_ dtt[Gr,{t)-Gr,{-t)], 

al = G'c(O), 

Thus, cr^ and can be explicitly determined from lattice calculations because the operators are evaluated at 
t = 0, which is located on the sampled path, <t <i/3. 

The remainder of this section is devoted to the prospects of determining r] and C- In particular, we investigate 
the usefulness of calculating moments of the path integral, 

CmnW)= / dhdt2V{Aih)A{t2))m. (48) 
Jo 
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For the remainder of the section, the subscripts a, 6, ij and C will be suppressed. Such moments can be calculated 
by considering a sampling of paths, each with a weight Wj. For each path, one can calculate "fj.m = J dt A(t)V" , 
which allows one to calculate 

c„„(ffl . ?a^^2i^. (4.,) 

The following analysis is aimed at evaluating how one might use Cmn to determine the viscosity. 
Using the fact that {A{ti)A{t2)) — G{t2 — ii), the moments can be expressed as: 

c„,r^{^3) = [''dtG{t)['' dhitTit+hr+t'^it+hr]. (50) 

Jo Jo 

The first several moments are: 

CooW = / dtGit)ii(3-t), (51) 
Jo 

CoiiPm = I dtG{t)i(3{i(3-t), 

ii3t{ip ~ t) + - tf 

o 

Thus, Coi provides no additional information. The reflection property of thermal propagators, 

G(t-i/3/2) = G(i/3/2-t), (52) 

means that G{t) has a reflection symmetry about the point t — i(3/2. Thus, only even moments of z = i — i/3/2 
are non-zero, and after substituting for t ~ z + i/3/2 in the expressions above, 

Coo(/3) = ^ dtG{t), (53) 

dt z^G{t), 



CiiiP) = / dtGit) 
Jo 

G20W) = / dtGit) 



and C20 provides no additional information about G{t). Thus, the moments Cmn appear sufficient to determine 
all the moments, 

DmiP) = -i dt z"G(<), z^t- i(3/2, (54) 
Jo 

for even m and Dm{(3) = for all odd m. The moments Dm(/3), which are real, contain all the information 
represented by the thermal propagators when calculated for imaginary times < i < Our goal is to 
determine the degree to which Dm{P) constrains the propagators for real time, and in particular, the integrated 
quantities in Eq. (|4T|) that represent the transport coefficients. 

As the transport coefficients involve integrating G{t) over real time, it is insightful to express the moments 
D„i in terms of integrals over real time. This can be accomplished by considering the dashed-line contour in 
Fig. [21 Whereas D„i represents the integral over the left-side boundary, integrating over the entirety of the 
boundary should give zero. Assuming SA is defined such that the propagator vanishes for large real times. Dm 
can be identified with the difference in the integrals over the upper and lower boundaries. Furthermore, the 
integral over the upper boundary can be identified with a integral over a path with negative real time using the 
reflection property of G{t) about the point t = i(3/2. Thus, 

/•oo 1^0 

Dm{(i) = -i I dt G{t){t-ip/2y'' + i{-iy'' dt G{t)it^if3/2y'' (55) 

^ — 00 

00 

dt [G{t){t - i(3/2r - G{-t){t + i(3/2r] 
dt {Gi{t)n{t - i/3/2)'" + Gfl(t)3(t - i/3/2)"} , 
Gj{t) ^ ^[G{t)-G{^t)], Gn^\[G{t) + G{-t)]. (56) 



11 




FIG. 2: (color online) Thermal path integrals can provide G{t) for t along the imaginary axis. Due to analyticity, the 
integral of G{t){t — ijS /2)™' along the contour (dashed line) is zero. Furthermore, the reflection properties, G(t — i/3/2) = 
G{iP/2 — t) (illustrated by the arrows), permit the integral along the imaginary axis from zero to if3 to be related to the 
difference of the integrals along the positive and negative real axes. 

The functions Gr and Gi correspond to the real and imaginary parts of G{t) for real t. It is instructive to view 
the last relation for the first few m, 

Do = 2 dt Gi, (57) 

JQ 

Di ^ 2 [ dt {tGi(t)^(p/2)GR{t)}, 
Jo 

D2 = 2 [ dt {{t^+l3''lA)Gi{t)-l3tGR{t)]. 
Jo 

From inspecting the expressions above, it is clear that the moments Di,D2 • • • , are insufRcient to completely 
determine all the moments of both G/ and Gr. In particular, the transport coefRcients, t] and C, are defined 
by the integral J dtGi{t)t, which is ambiguous given that Di also depends on Gr. The inability to determine 
all the moments of interest is expected given that G{t) was only determined for < t < if3. If G had been 
known for all t along the imaginary axis, it would permit the determination of the Fourier transform of G, 
which could then be used to determine G along the real axis by analytic continuation. In Reference [4^] the 
expression for Dq in Eq. (|57|) is exploited to make a link between lattice results and viscosity, by assuming that 
/p dt Gi is related to the viscosity by a characteristic relaxation time, assumed to be the inverse temperature. 
However, near Tc the relaxation time could be much larger, since there is little thermodynamic motivation for 
the system to move towards the absolute equilibrium given the broad minimum to the free energy associated 
with competing phases [sj. 

Given the limitations discussed above, determining the transport coefficients will rely on assuming a functional 
form for G{t), whose parameters can be fit to the moments DmiP)- The functional form can then be applied to 
determine the transport coefficients. In the next three sub-sections three forms are explored. The ffist is based 
on Gaussians modified by Hermitc polynomials, and the second on an exponential form. The third form [TsLflGj 
is based on a Breit-Wigner form motivated by a single mode which oscillates with a given frequency and width. 
An alternative to using moments of G{t) to determine parameters for a fitting function is to use G(t) at specific 
values of t along the imaginary axis. These are also addressable with path integrals by setting t to the discrete 
values represented by the lattice. 



A. Expanding G{t) with a Gaussian and Hermite polynomials 

Given the symmetry constraints, G{t — j/3/2) = G{—{t — i/3/2)), and the constraint that G must vanish for 
large real time, Gaussians, and the associated orthogonal polynomials seem a reasonable basis to expand, 

G(t) =e-^'/^' [go + g2H2{z/T) + gMz/r) ■ ■ ■]e-''"^^^\ z = t-i(3/2, (58) 

where Hi are Hermite polynomials, with the odd polynomials being set to zero from symmetry. The first few 
even polynomials are: 

Ha = 1, Hi{x) = 2x, H2{x) = Ax^ - 2, H:i{x) = 8a;^ - 12a;, Hi{x) ^ IQx'^ ~ A^x^ + 12, (59) 
X = \hi{x), x^^\ {H2(x) + 2Ho{x)) , ^3 = i (H:i{x) + 6ili(x)) , a;^ = (^4(a:) + \2H2ix) + \2Ha{x)) . 

2 4: o io 
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Given the polynomials, one can perform the integrals in Eq. (|54p to determine the momentsDm (/?) , 



* / dz V 5„z'"7J„ (z) exp f-L- , 



(60) 



By expanding z^Hm in powers of z using the first line of Eq. ([59|) . then rewriting each term as a sum over Hm 
using the second line of Eq. (|59p . one then can perform the integrals by exploiting the relations, 

(61) 



In{x) = j dy e y Hniy), 

In>o{x) = -i Hn-i{Q) ~ e^^ Hn-i{ix) 
If the expansion of G{t) is second order, i.e. gn>i = 0, 
-Do = T(.go/o + 32/2)6"'=', 

D2 = -y {go (/o + /2/2)+. 92 (4/0 + 5/2+14/2)} e-"', 



(62) 



Da 



— {go (3/0/2 + 3/2/2 + 74/8) + .92 (9/0 + 33/2/2 + 7/4/2 + Ie/16)} 



where /„ are evaluated at /„(x = /3/2t). Thus, if lattice calculations were to extract Dq, D2 and D4, one could 
then find the three parameters r, 170 and g2. Including the g^ term in the expansion of G{t) would then require 
Dq to constrain the parameters. 

The simplest expansion would be to keep only the zero'th order Hermite polynomial, i.e., gn>2 = 0. In that 
case, one would only require knowledge of Dq and D2 to determine the two parameters go and r. For this case, 
as can be seen from inspecting Eq. (j62p . x = /3/2t is determined by the ratio, 



{pm^Do 



2Io{x)+l2{x) 
x'^Ioix) 



(63) 



This allows x to be extracted by inverting the ratio. However, the r.h.s. of the equation never falls below 1/3 
if only zero*^ order polynomials are used. Thus, if lattice results were to yield a lower ratio for the left-hand 
side, the expansion of G(t) would require more Hermite polynomials. This limit comes from the fact that, when 
keeping on the zero*'' term, G(t) has a minimum at < = i/3/2 when traveling along the imaginary axis. This 
constraint is lifted by keeping additional terms. 

The transport coefficients are defined by the integrals along the real axis of G/, the imaginary part of G{t), 
as defined in Eq. (|5S|) . The viscosities become (using 77 to denote either the shear or bulk viscosity) 



- dt tGi{t) - - / dt tG{t) 
^ J -00 ^ J —00 



(64) 



oo-i/3/2 



dz (z + il3/2) [goHo{z/T) + .92//2(z/t) + giHiiz/r) ■ 



-(z=^+/3V4)/t 



where z = t — i(3/2. From analyticity, the integral will not be changed by altering the limits on the path to ±00. 
After recognizing (z + i/3/2) — [tHi +i(3Ho)/2, one can use orthogonality properties of the Hermite polynomials 
then allow one to ignore all the terms except those involving Hq and obtain 



The fluctuations become 



a' = G(0) = 50 + 52/^2 (-1/3/2T) + giHi{^il3/2T) + 



(65) 
(66) 



The fluctuation G{t — 0) can also be calculated directly from lattice calculations, without taking moments. 

One can see from Eq.s ((65)l and ([66)1 that the ratio r^/cr^ is limited if one were to expand Go with only zero*'' 
order Hermite polynomials. This is related to the inability of such an expansion to cover all possibilities of 
D2/ Do as mentioned earlier. 

As with most expansions, the validity is determined by the convergence. If the resulting transport coefficients 
do not appreciably change as the number of Hermite polynomials are increased, the results should be considered 
robust. 
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B. Exponential expansions 

The most intuitive form for G(t) would involve exponential decay for large times. Such a form was considered 
in [il. Choosing such a form to satisfy the reflection properties around t = i(3/2, 

AT f , cosuj(t-iP/2) 

The result has the desired exponential behavior, but is non-analytic along the imaginary axis, precisely where 
the thermal path integrals are evaluated. One can forge ahead and calculate the moments Dm by assuming that 
Git) = (1/2)(G(< + e) + Git - e)) for m = 0, 

/•/3/2 

D,n=0 2i-^A dx iixy^cosiVx). (68) 

J-P/2 

Due to the non-analyticity of G{t) the odd moments, which are zero when integrating along the imaginary axis, 
are non-zero when integrating an amount e on either side of the axis. However, the calculation of the even 
moments are unchanged by the infinitesimal translation. Thus, only the m = 0,2,4 - • • thermal moments can 
be used to constrain the corresponding moments along the real axis. 
The first two even-numbered moments are: 

2A AA 
i?o = — sin(/3r/2), i?2 = -_sin(/3r/2). (69) 

The decay width T is thus obtained by taking the ratio of the first two moments, 

r^ = -f^. (70) 



The viscosities are straight-forward to find by integrating the form of Eq. ([GT]) over real time as was done in 

Eq. dell), 

2A 

r?=— sin(/3r/2). (71) 

One could increase the complexity of the form for Git) in Eq. (|67p by incorporating an expansion in expo- 
nentials with width F, 2r • • • , with a similar tact as was used for Hermite polynomials in the preceding section. 
However, it is difficult to assess the inherent error in using a form which is expressly non-analytic. It may well 
turn out that the exponential form is valid if the characteristic lifetime, l/F is much larger than /?. This corre- 
sponds to the classical limit, in that Vfi ^ as ?i ^ 0. Furthermore, if one calculates the effective lifetime, r^/a^ 
for the Gaussian/Hermite expansion in Eq.s (j65p and (j66p . it appears that that expansion might have difficulty 
in the opposite limit. Thus, different forms might become preferable depending on whether the characteristic 
lifetimes are large compared to the inverse temperature. 



C. Breit-Wigner Form 

In [l^ the shear viscosity was determined from lattice calculations by investigating G(t) for imaginary t and 
fitting to a Breit-Wigner form proposed in [l^. This form has a physical motivation in that G(t) corresponds 
to a thermal Green's function for a single mode. Using simple Bose creation/destruction operators. 

Git) = A{ia + o))iait)^a\t))) (72) 
= A [/(m)e-""* + (1 + /(m))e""*] , 
- A [/(m)e-™* - /(-m)e™*] . 

where /(m) = (a^a) = 1/(6^"* -f 1) is the Bose occupation factor. This can be rewritten in terms of a spectral 
function p(cij) = ASiuj — m) — A5i-~L0 — m), 

Git) = fdoj -4-^p(c.)e-^-*. (73) 
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The Karsch-Wyld assumption [1^ for p(w) is a Breit-Wigner form, 

7 7 



TT 



(to — oj)^ + 7^ (m + a;)^ + 7^ 
It is straight-forward to calculate the viscosity from this form. 



(74) 



Ti = hm ——^ i (75) 



(7^ -f- TO^)2 

In order to calculate the moments Dm-, one can rearrange G{t) to the form, 

^/ N /" , cos( Lult — iB/2)) , , 
7 smh(/3Ci;/2) 

By comparing this form to the equivalent form for the exponential form in Eq. (|67p . one can see that the two 
expressions differ by the substitution 

1 , P(-) (77) 



cj2 + r2 sinh(/3cj/2) 

To finish finding D„n one inserts the expression for G{t) into the definition for D„i in Eq. (|54p and obtain 

/ , P{^) f^'^ , w X 2Ato 



f , P 0^ [ ' , w N 2Ato 

/ . , )J / cosh = ^ ^, (78 

J smh /5t^/2 7_«/2 ' ' to2 -h 72 ' ^ ' 



D2 = — I duj , , f dz cosh(w2;) 

' sinh(/3c./2) 1^/2 ^ ^ 



/ ■ ■ ._9Cosh(/3(m -I- 17)72) 1 

smh(/5(TO -I- i7)/2) J 



^ 87rAp(a;„) 

-1-^ ^ ^ , i^n — In-Klj p. 

n— 1,00 

Since the Karsch-Wyld form is based on three parameters (A, m and 7), one must either calculate D^, or use 
ct2 = G(0). 

CT^ = 2A5R{coth(/3(TO + i7)/2)}- ^3 ^ p(t^„)- (79) 

Tl — 1,00 

Rather than fitting the moments, and G(t = 0), one could alternatively find the best fit to G(t) evaluated at 
the lattice points, as was done in p^ . 

Like the exponential form of the previous section, the Breit-Wigner form is non-analytic. This can be seen 
by calculating Gi{t), which unlike Gr, has a simple form. From Eq. ([73]), 

Gi{t) = Y (^(^) - ^(-*)) ^ ^ / P^^'^ 

= 2Ae^^* sm{mt)Q{m) + 2Ae^* sm{mt)e{-m). 

Although the form is not analytic for t along the imaginary axis, it is continuous, unlike the exponential form 
of the previous section. This maintains the condition -Dm=odd = for the path illustrated in Fig. [21 The most 
distinct difference that separates the Breit-Wigner from the exponential form of the previous section is in the 
oscillation described by the sin(TOt) term in Eq. (|80p . An example where such an oscillation might be physical 
concerns the bulk viscosity arising from a non-equilibrium mean field, which might be under-damped. The trace 
of the stress-energy tensor might then oscillate along with the field. It is more difficult to motivate oscillatory 
behavior for the shear terms in the stress energy tensor (the terms in the first section), which in the classical 
limit, should decay exponentially. One way to mimic exponential behavior would be to set to = 0, which would 
lead to Gi{t) ~ te-^K 



15 





Hermite 


Exponential 


Breit-Wigner 


G{z) 
parameters 
Do 
D2 

V 


T, 90,92 ■ ■ ■ 
2re-^'E„5n/„W, x = l3/2r 
Eq. (inij 
E„ 9nH^{-ix) 
(0F/2)/3r<7oe-^' 


A,r 

2^sin(/3r/2)/r 
-4Asin(/3r/2)/r^ 

^cos(/3r/2) 
2Asin(/3r/2)/r2 


A / dujp{ij)e-''^y sinh(/?tj/2) 

2y4m/(m^ +7^) 
Eq. JZHl) 
Eq. ([79)1 

4Am7/(7^ + m^)^ 



TABLE I: Three functional forms for G(z = t — j/3/2) are given above. In principle, lattice calculations can explicitly 
provide the variances for the shear- and bulk-related components of the stress-energy tensor, as well as moments of 
the Green's function, Dm defined in Eq. (I54|l , integrated along the imaginary axis. Fitting either these moments, or the 
actual values of G at the lattice points, to a given functional form, can then provide a value for the viscosity. However, 
the value could vary substantially from form to another. 



D. Summary of Parameterizations 

Path integrals performed along the imaginary axis can explicitly provide and cr^, which when given the 
viscosities 17 and ^, yield the IS relaxation times through Eq. (jl5p . Determining 77 and ^ requires fitting to a 
functional form, the choice of which should be motivated by the expected physical behavior. For instance, if 
one does not expect oscillatory behavior of the relevant elements of the stress-energy tensor, one should avoid 
the Breit-Wigner form. A second consideration concerns the analyticity of G(t) for t along the imaginary axis. 
Since the time dependence is driven by the evolution operator, e*^*, it should be analytic. Furthermore, the 
characteristic times for any changes in G{t) should not be shorter than microscopic time scales. It is difficult to 
assess the dangers of using explicitly non-analytic forms, such as the exponential and Breit-Wigner forms, but 
intuitively one would expect the errors to be important if the time scales are small. For long time scales, one 
would not expect non-analyticities, which arc unphysical if one is considering the short-time behavior of G/(t), 
to be important. 

One limit in which the various parameterizations can be compared is the classical limit. In that limit (l^ . 

e^^" Ait)e-'^" « A{t) + if3hdtA{t) + Oh^, (81) 
which allows one to approximate Gj as 

- i{[A,A{t)]) ^ f3dt{AA{t)) + Oh. (82) 

With this approximation, 

V = dt t{[A, A{t)]) K (3 dtG{t), (83) 

Jo Jo 

Do = dt {[A,A{t)])^/3a\ 
Jo 

The classical approximation works whenever f3 is sufficiently small that when multiplied by characteristic mi- 
croscopic energy scales, one should get a number much less than unity. This should be true for either large 
temperatures or long relaxation times, i.e., Tt » h. The former expression is satisfied by all three forms in 
Table H] as /3 — > 0. One way to check the validity of the classical limit is to compare to Dq, both of which are 
calculable on the lattice. 

In the classical limit, all three forms provide different answers for the viscosities given the same values for 
Dm- For instance, if one uses a strictly Gaussian expression for the Hermite expansion, gm>o = 0, one can find 
that as 77 ^ 0, 77 ^ Z3o-D27r/24, whereas for the exponential form 77 ^/—DoD2/2. Thus given the same 
values Do and D2 from the lattice, the two forms would lead to viscosities that differ by nearly a factor of two. 
The Breit-Wigner form can vary even more, since the addition of a third parameter makes 77 under-constrained 
if one uses only Do and D2- Physically, the oscillatory nature of the Breit-Wigner form makes the moments less 
predictable, and underscores the importance of selecting a non-oscillatory functional form if there is physical 
justification. Of course, the Hermite expansion can mimic oscillatory behavior of G{t) if given a sufficient 
number of terms. 

One feature of all three functional forms described above is that they can accommodate arbitrarily short 
relaxation times. In some of the arguments used to motivate the KSS limit, 77/s > Tl/Att [iTllTsj. appeals have 
been made to the energy/time uncertainty principle in the form Tt ^ h. This constraint is easily violated by 
all three functional forms considered here. 
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VII. SUMMARY 



The considerations presented here all seem to solidify Israel-Stewart approaches as the preferred approach for 
hydrodynamic descriptions of high-energy heavy-ion collisions. Israel-Stewart treatments are more numerically 
stable and physical. Although the treatments require additional parameters, they are uniquely determined by the 
microscopic properties of the medium, which are formally stated in a way that can be addressed in sophisticated 
approaches like lattice gauge theory. Israel-Stewart parameters are only uniquely determined for small deviations 
from equilibrium, even if the velocity gradients are large, but treatments can be phenomenologically extended 
to describe large deviations by adding extra parameters to confine the magnitude of the deviation of the stress- 
energy tensor to physical values. 

Since lattice calculations do not directly explore thermal expectations of operators evaluated at real times, 
neither the viscosities or relaxation times are directly accessible. However, the ratios of the viscous parameters 
to the relaxation times, which correspond to equal-time thermal fluctuations of the stress-energy tensor, can 
be explicitly determined. Extracting the viscosities requires fitting the behavior of the thermal correlators, 
which can be directly determined for imaginary times, < t < i(3, to a chosen functional form, which is then 
analytically continued to real times. The only firm lessons to be learned from the comparison of the three 
functional forms investigated in Sec. IVII is that the choice of form can significantly affect the answer, and 
that the choice of form should depend on whether one expects oscillatory behavior, and whether the classical 
approximation (long relaxation times compared to the inverse temperature) is justified. 
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